%create synthetic mogi waveforms

fs = 100;
tspan = 60*60;
time = (0:1/fs:tspan).';
if mod(length(time),2) == 1
    time = time(1:end-1);
end
f = (fs*(0:ceil(length(time)/2))/(length(time))).'; %interpolated one-sided freq (HZ)

 %%%%import seismic data
    %import HVO station locations
    %stations = {'KKO'; 'BYL'; 'HAT'; 'SBL'; 'UWB'; 'NPT'; 'OBL'; 'UWE'; 'SDH'; 'WRM'; 'RIMD'; 'NPB'; 'SRM'};
    stations = {'KKO'; 'BYL'; 'HAT'; 'SBL'; 'UWB'; 'NPT'; 'OBL'; 'UWE'; 'SDH'; 'WRM'; 'RIMD'};
    stationindex = true(size(stations));
    stations = stations(stationindex);
    locations = cell(size(stationindex));
    locations{1} = struct('station','KKO', 'latlon', [19.39781 -155.26600 1146.0]);
    locations{2} = struct('station','BYL', 'latlon', [19.41200 -155.26000 1059.0]);
    locations{3} = struct('station','HAT', 'latlon', [19.42300 -155.26100 1082.0]);
    locations{4} = struct('station','SBL', 'latlon', [19.42700 -155.26800 1084.0]);
    locations{5} = struct('station','UWB', 'latlon', [19.42469 -155.27800 1091.0]);
    locations{6} = struct('station','NPT', 'latlon', [19.41200 -155.28100 1115.0]);
    locations{7} = struct('station','OBL', 'latlon', [19.41750 -155.28400 1107.0]);
    locations{8} = struct('station','UWE', 'latlon', [19.42111 -155.29300 1240.0]);
    locations{9} = struct('station','SDH', 'latlon', [19.38950 -155.29400 1133.0]);
    locations{10} = struct('station','WRM', 'latlon', [19.40650 -155.30000 1163.0]);
    locations{11} = struct('station','RIMD', 'latlon', [19.39531 -155.27400 1128.0]);
    %convert lat-lon locations to UTM E-N coordinates
    zone = utmzone(locations{1}.latlon(1),locations{1}.latlon(2));
    utmstruct = defaultm('utm');
    utmstruct.zone = zone;
    utmstruct.geoid = wgs84Ellipsoid;
    utmstruct = defaultm(utmstruct);
    for i = 1:length(locations)
        locations{i}.loc = [0 0 0];
       [locations{i}.loc(1), locations{i}.loc(2), locations{i}.loc(3)] = mfwdtran(utmstruct,locations{i}.latlon(1),locations{i}.latlon(2),locations{i}.latlon(3)); 
    end
    
    %import waveforms and instrument responses
    s = cell(size(stations));
    count = 0;
    for i = 1:length(locations)
        if ismember(locations{i}.station,stations)
            count = count + 1;
            s{count} = struct('station',locations{i}.station);
            s{count}.time = time;
            if mod(length(s{count}.time),2) == 1
                s{count}.time = s{count}.time(1:end-1);
            end
            s{count}.loc_east = locations{i}.loc(1);
            s{count}.loc_north = locations{i}.loc(2);
            s{count}.loc_vert = locations{i}.loc(3);
            %instrument responses
            data = csvread(strcat(s{count}.station,'_vel_response.txt'),1);
            s{count}.raw_response = data(:,2)+1i*data(:,3);
            s{count}.raw_f_response = data(:,1);
            %interpolate instrument response
            s{i}.response = interp1(s{i}.raw_f_response,real(s{i}.raw_response),f) + 1i*interp1(s{i}.raw_f_response,imag(s{i}.raw_response),f);
            s{i}.f_response = f;
            response(i,:) = s{i}.response;
            response(length(s)+i,:) = s{i}.response;
            response(2*length(s)+i,:) = s{i}.response;
        end
    end
    %
    
    %%%%lava lake location (from Chouet work)
    epilatlon = [19.405402 -155.281041 100.0];
    [epiloc_east, epiloc_north, epiloc_vert] = mfwdtran(utmstruct,epilatlon(1),epilatlon(2),epilatlon(3));
    %}

    
    %%%% parameters
    mu = 30*10^9; %shear modulus
    nu = 0.25; %poisson's ratio
    %inversion parameters in form [min,max,guess]
    eastoffset = 0; %centroid offset from epiloc if dip is positive
    northoffset = 0; %centroid offset from epiloc if dip is positive
    depth = 600; %should be positive 
    
    centroideast = eastoffset+epiloc_east;
    centroidnorth = northoffset+epiloc_north;
    
    east = centroideast;
    north = centroidnorth;
    
    obs = zeros(3,length(s)); %observation (station) locations
    obstilt = zeros(3,4*length(s));%observation locations for tilt
    for i = 1:length(s)
        obs(1,i) = s{i}.loc_east;
        obs(2,i) = s{i}.loc_north;
        obs(3,i) = 0;
        obstilt(:,4*(i-1)+1) = obs(:,i) + [-1;0;0];
        obstilt(:,4*(i-1)+2) = obs(:,i) + [1;0;0];
        obstilt(:,4*(i-1)+3) = obs(:,i) + [0;-1;0];
        obstilt(:,4*(i-1)+4) = obs(:,i) + [0;1;0];
    end 
    

    for j = 1:length(s)
        s{j}.east_pred = zeros(size(s{j}.time));
        s{j}.north_pred = zeros(size(s{j}.time));
        s{j}.vert_pred = zeros(size(s{j}.time));
    end

    for j = 1:length(s)
        s{j}.east_pred_tilt = zeros(size(s{j}.time));
        s{j}.north_pred_tilt = zeros(size(s{j}.time));
        s{j}.vert_pred_tilt = zeros(size(s{j}.time));
    end

    

        %calculate greens functions for inflate-deflate motion
        mdl = [depth east north 100]';
            U = mogi_greensfunc(mdl,obs,mu,nu);
            Utilt = mogi_greensfunc(mdl,obstilt,mu,nu);
        G = [U(1,:).';U(2,:).';U(3,:).'];
        Gtilt = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(s),1)];
%     end

%synthetic up force
% %up = sin(2*pi*time(1:end-2000)/30);% + sin(2*pi*time(1:end-2000)/20);
% %up = up.*sin(2*pi*time(1:end-2000)/(time(end-2000)))*10^2;
% % up = 0*time(1:end-2000);
% % up(floor(end/4):floor(3*end/4)) = 1;
% up = (0:length(time(1:end-2000))-1).';
% sigtaper = 10^-2;
% up = 1./(1+exp(-sigtaper*(up-length(up)/4))).*1./(1+exp(sigtaper*(up-3*length(up)/4)));
% up = [zeros(1000,1);up;zeros(1000,1)];

up = 0*time;

T = 40;
Q = 25;
A = 2E5;
offset = 10*60;
uptemp = A*sin(time*2*pi/T).*exp(-time*2*pi/Q/T);
uptemp = [zeros(offset*fs,1);uptemp(1:end-offset*fs)];
up = up+uptemp;

T = 15;
Q = 40;
A = 1E5;
offset = 9.8*60;
uptemp = A*sin(time*2*pi/T).*exp(-time*2*pi/Q/T);
uptemp = [zeros(offset*fs,1);uptemp(1:end-offset*fs)];
up = up+uptemp;

T = 12;
Q = 45;
A = 1E5;
offset = 9.81*60;
uptemp = A*sin(time*2*pi/T).*exp(-time*2*pi/Q/T);
uptemp = [zeros(offset*fs,1);uptemp(1:end-offset*fs)];
up = up+uptemp;

% Anoise = 1E4;
% up = up + Anoise*randn(size(up));
% !!adding noise below!!

Fop = fft(up);
Fop = Fop(1:end/2+1);
Fop(end) = conj(Fop(end));


    %use greens function matrix to calculate best fit [ss,ds,op] at each time
    Fss = zeros(length(f),1);
    Fds = zeros(length(f),1);
    %Fop = zeros(length(f),1);
    FU_pred = zeros(3*length(s),length(f));
    FU_pred_tilt = zeros(size(FU_pred));
    g = 9.807; %leave positive
    for i = 2:length(f)
%         if f(i) < 4*f_lowpass %skip high frequencies that have been filtered out anyway
%             if onlyop
%                responsef = response(:,i);
%             else
%                responsef = [response(:,i),response(:,i),response(:,i)];
%             end

            %M = (G)\FU_real(:,i); %no tilt
            %FU_pred(:,i) = (G)*M; %no tilt

            %ir removed earlier
            %M = (G + g*Gtilt/(2*pi*1i*f(i))^2)\FU_real(:,i);
            M = [Fop(i)];
            FU_pred(:,i) = (G + g*Gtilt/(2*pi*1i*f(i))^2)*M;
            FU_pred_tilt(:,i) = (g*Gtilt/(2*pi*1i*f(i))^2)*M;

%             M = (responsef.*(G + g*Gtilt/(2*pi*1i*f(i))^2))\FU_real(:,i);
%             FU_pred(:,i) = (responsef.*(G + g*Gtilt/(2*pi*1i*f(i))^2))*M;
%             FU_pred_tilt(:,i) = (responsef.*(g*Gtilt/(2*pi*1i*f(i))^2))*M;

%             if onlyop
%                 Fop(i) = M;
%             else
%                 Fss(i) = M(1);
%                 Fds(i) = M(2);
%                 Fop(i) = M(3);
%             end
    end
    
    %Multiply with response
    fmat = repmat(f.',33,1);
    %!!!disp_resp = vel_resp.*(1i*2*pi*f), seems odd but otherwise IR
    %removal doesn't make sense (ie, counts/vel_resp = vel and
    %counts/(vel_resp*f) = disp
    FU_pred = FU_pred.*response.*(1i*2*pi*fmat);
    FU_pred(:,1) = 0;
    
    %make 2-sided spectra
    FU_pred = [FU_pred(:,1:end-1),fliplr(conj(FU_pred(:,2:end)))];
    
    %Calculate U_pred
    U_pred = zeros(3*length(s),length(s{1}.time));
    for i = 1:3*length(s)
        U_pred(i,:) = real(ifft(FU_pred(i,:)));
        %U_pred(i,:) = U_pred_temp(end/3+1:2*end/3);
    end
    
    %add noise
    U_pred = U_pred + 0.01*range(U_pred(:))*randn(size(U_pred));
    
    for i = 1:length(s{1}.time)
        for j = 1:length(s)
            s{j}.east_pred(i) = U_pred(j,i);
            s{j}.north_pred(i) = U_pred(length(s)+j,i);
            s{j}.vert_pred(i) = U_pred(2*length(s)+j,i);
        end
    end
    
    'not saving'
    %save output
%     for i = 1:length(s)
%         matrix = [s{i}.time,s{i}.east_pred,s{i}.north_pred,s{i}.vert_pred];
%         filename = strcat('halemaumau_mogi_step_synthetic_',s{i}.station,'.txt');
%         dlmwrite(filename,matrix,'precision',10)
%     end
    
    